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Abstract 



In the functional approach to quantum chromo dynamics, the properties of hadronic bound states are accessible via 
' covariant integral equations, e.g. the Bethe-Salpeter equations for mesons. In particular, one has to deal with linear, 
homogeneous integral equations which, in sophisticated model setups, use numerical representations of the solutions of 
other integral equations as part of their input. Analogously, inhomogeneous equations can be constructed to obtain off- 

C" l shell information in addition to bound-state masses and other properties obtained from the covariant analogue to a wave 
function of the bound state. These can be solved very efficiently using well-known matrix algorithms for eigenvalues 
|(in the homogeneous case) and the solution of linear systems (in the inhomogeneous case). We demonstrate this by 
solving the homogeneous and inhomogeneous Bethe-Salpeter equations and find, e.g. that for the calculation of the mass 

qq spectrum it is more efficient to use the inhomogeneous equation. This is valuable insight, in particular for the study of 
baryons in a three-quark setup and more involved systems. 
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The underlying quantum field theory of the strong inter- 
action in the standard model of elementary particle physics 
is quantum chromodynamics (QCD), a non-abelian gauge 
theory which deals with elementary degrees of freedom 
called quarks and gluons [lj]. A remarkable feature of QCD 
is asymptotic freedom, which means that the runningcou- 
pling of the theory is small in the high-energy regime [2H4| . 
There, perturbation theory can be applied, and pertur- 
bative QCD has been well established in the high-energy 
domain, (e.g. [B| and references therein). At low ener- 
gies, however, perturbation theory is no longer applicable, 
since the value of the running coupling increases to the 
order of 1. Since bound states are intrinsically nonpertur- 
bative, corresponding methods have been developed and 
used to investigate hadrons, the bound states of quarks 
and gluons. We eclectically list a few references regard- 
ing constituent quark models foj- TToll . effective field theories 
, lattice-regularized QCD 14- III, QCD sum rules 
, and the renormalization-group approach to QCD 
27j (always see also references therein). 
Another remarkable property closely related to bound 
states is the so-called confinement of quarks and gluons. 
It entails that only objects like hadrons, where the color 
charges carried by the elementary degrees of freedom are 
combined to a color-neutral state, can be observed directly. 
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While in constituent-quark models confinement is usually 
implemented via potential terms of an infinitely rising na- 
ture (of, e.g., harmonic-oscillator or linear type), in QCD 
the particularities are more delicate (for a recent review of 
the problems surrounding quark confinement, see e.g. [28}). 
In a quantum field theoretical setup, as we use it here, con- 
finement is tied to the properties of the fundamental Green 
functions of the theory. 

In the present work, we employ the Dyson-Schwinger- 
equation (DSE) approach to QCD. The DSEs are the co- 
variant and nonperturbative continuum equations of mo- 
tion in quantum field theory. They constitute an infinite 
set of coupled and in general nonlinear integral equations 
for the Green functions of the quantum field theory un- 
der consideration. There are several extensive reviews on 
the subject that focus on different aspects of DSEs, like 
fundamental Green functions [29 31], bound-state calcu- 
lations [H, |33[ and applications of the formalism, e.g. to 
QCD at finite temperature and density [34| . Bound states 
are studied in this approach with the help of covariant 
equations embedded in the system of DSEs. In particu- 
lar, the Bethe-Salpeter equation (BSE) ( 35 [. l36| is used for 
two-body problems such as mesons [37H39l | and covariant 
Faddeev-type equatio ns 1401 are used for three-body prob- 
lems such as baryons [4l| . 

Ideally, one could obtain a self-consistent simultaneous 
solution of all DSEs, which would be equivalent to a so- 
lution of the underlying quantum field theory. While in 
investigations of certain aspects of the theory such an 
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approach is successful (see, e. g. [42|, |43j and references 



therein) , numerical studies of hadrons necessitate a trun- 
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cation of this infinite tower of equations. 

Once the covariant bound-state equation has been 
solved to obtain the mass spectrum of the system un- 
der consideration, the corresponding covariant amplitudes 
can be used to compute further observables. In rainbow- 
ladder truncation, prominent examples include leptonic 
decay constants [43446} , hadronic decays [ 47l[48l , and elec- 
tromagnetic properties of both mesons [49l453j and baryons 



geneous (vertex) BSE, 

r [2] (fc,p) = r (fc,p) 
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59j . Improvements to this truncation have been con- 
sidered in the past and studies in this direction are un- 
der way [60l465l | . What we discuss in the present work 
is most easily exemplified in a simple truncation, but be- 
comes more important — and thus relevant — with any 
kind of increasing numerical effort necessitated by either a 
more involved truncation or the study of a system of more 
than two constituents. 

The paper is organized as follows: in sec. [5] we collect 
the necessary formulae regarding covariant bound-state 
equations as they are obtained in the DSE approach to 
QCD. Section |3] details the discretization of the integrals 
and the general numerical setup. Section [4] contains nu- 
merical solution strategies for both the homogeneous and 
inhomogeneous bound state equations. In sec. [5] we ap- 
ply the methods described to solve the homogeneous and 
inhomogeneous BSE for pseudoscalar mesons in rainbow- 
ladder truncation and analyze the efficiency of the algo- 
rithms. Conclusions and an outlook indicating both im- 
mediate and further possible applications of the strategies 
described herein are offered in sec. [SJ 

All calculations are performed in Euclidean space. 



2. Structure of covariant bound-state equations 



In the DSE approach to QCD, mesons are described by 
general vertices connecting (anti-)quarks to objects carry- 
ing the appropriate quantum numbers as demanded by 
the respective superselection rules. These vertices are 
the so-called (inhomogeneous) Bethe-Salpeter amplitudes 
(BSAs), denoted by Tt 2 ](k,P), which describes a two- 
particle system, denoted by the subscript pi, with to- 
tal momentum P and relative momentum k of the con- 
stituents. The inhomogeneous BSA satisfies the inhomo- 



+ K [2] (k,q,P)S a (q + )T m (q,P)S b {q-) , (I) 

Jq 

where T°(k,P) is a driving term with the quantum num- 
bers of the system, the Euclidean-space four-dimensional 
momentum integration is given by J q = J tj^t , S a ' b (q±) 
denote the renormalized dressed (anti-)quark propagators, 
K[2] {k, q, P) represents the quark-antiquark interaction 
kernel, and q± — q ± rj±P are the (anti-)quark momenta 
with momentum partitioning r/± such that t] + +t]- = I, re- 
spectively. The choice of the momentum partitioning is ar- 
bitrary and usually a matter of convenience, e.g., for equal- 
mass constituents a convenient choice is i] + = r\- = 1/2. A 
graphical representation of Eq. (fTJ) is given in Fig. [TJ where 
the arrows denote dressed-quark propagators (analogously 
in Figs. [5] and El). 

The solution of ([T]), T^(P,k), contains both off-shell 
and on-shell information about the states in a channel 
with the quantum numbers under consideration, which 
are fixed via the construction of T°(k,P) and T[ 2 ](P, k). 
In particular, T[ 2 ]{P, k) has polesJII whenever the total- 
momentum squared corresponds to the square of a bound- 
state mass in this channel (e.g. [(36[ and references therein). 
If there exists a bound state and the corresponding on- 
shell condition, in Euclidean space P 2 — —M 2 , is met, 
the properties of the bound state are described by the 
pole residues of Eq. (fT]). These residues, the homogeneous 
BSAs T[2h](k, P), can be obtained from the corresponding 
homogeneous BSE, 



r [2h ](k,P) = / K l2] (k,q,P)S a (q + )r [2h] ( q ,P)S b (q-) , 
Jq 

(2) 

depicted in Fig. [2] 

For Baryons, an analogous construction can be made, 
where the homogeneous equation for the on-shell ampli- 
tude is a covariant three-quark equation often referred to 



1 lt should be noted that this formalism is also applicable to res- 
onances, where instead of a pole in rpi (P, k) as a function of a real 
variable P 2 one expects a finite peak. 
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Figure 2: The homogeneous BSE, Eq. J3J 



Figure 3: The homogeneous equation for a three-body bound state 
(covariant Faddeev equation), Eq. J3]l 



as a covariant Faddeev equation 4l|, |67 
be written as 



, which may 



T l3h] (h,k 2 ,P)= I S a ( Pl )S b ( P2 )S c ( P3 ) 

x AT[ 3 ](fci, k 2 , qi, (? 2 , P)F [3h] (qi, q 2 , P) , (3) 

and a pictorial representation is given in Fig. [3J Here, 
the kernel K\s\(ki, k 2 , q\, q 2 , P) subsumes all interactions 
of the three quarks with the individual momenta pi, i — 
1,2,3, and the bound state is described by the covariant 
three-quark on-shell amplitude Tt 3 }^(ki,k 2 , P), which de- 
pends on the total momentum P as well as two relative 
( Jacobi) momenta k\ and k 2 . Note that this equation con- 
tains an integral over two momenta, namely q\ and q 2 , thus 
inflating the size of the problem in terms of a numerical 
setup. 

While in this work we will focus on the solution of 
bound-state equations such as (JlJ or ©, a note on the 
construction of the interaction kernels K and the origin 
of the quark propagators S as inputs in these equations 
is in order. In QCD, in addition to the set of DSEs, the 
Green functions of the theory also satisfy Ward-Takahashi- 
and/or Slavnov- Taylor identities, e.g. |29l. |30|. These re- 
late certain Green functions among each other and pro- 
vide guidance or even constraints in many cases, if one 
is to use a truncation and wants to make an Ansatz for, 
say, Km\. For light-hadron physics, the axial- vector Ward- 
Takahashi identity is of particular interest, since it en- 
codes the chiral symmetry of QCD with massless quarks 
as well as its dynamical breaking (see, e.g. [44], [63] for 
details). In other words, satisfaction of this identity guar- 
antees that the properties of the pion, the lightest hadron 
and would-be Goldstone boson of dynamical chiral sym- 
metry breaking follow the expected pattern, e.g. the pion 
mass vanishes in the chiral limit, and leads to a general- 
ized Gell-Mann-Oakes-Renner relation valid for all pseu- 
doscalar mesons 6^, 7(|. The rainbow-ladder truncation 
of the DSE-BSE system satisfies the axial-vector Ward- 
Takahashi identity. Beyond rainbow-ladder truncation, 
satisfaction of this identity can be achieved on more gen- 
eral terms and has been implemented throughout light- 
hadron studies of the past years in this approach, e.g. [60|— 

H, [zH-Izi] . 

More concretely, the satisfaction of this identity leads 



to a close relation of the interaction kernel Ki 2 ] in the 
quark-antiquark BSE and the quark self-energy present in 
the DSE of the quark propagator, the QCD gap equation. 
The consistent use of such related kernels in the BSE and 
the gap equation provides the proper input for the BSE 
in terms of the quark propagators S as solutions of the 
symmetry-preserving version of the gap equation. This 
consistency can be maintained also in numerical studies 
with high accuracy. In the following we always assume 
that the gap equation has been solved with the appropri- 
ate self-energy to match K[ 2 ] and that the resulting quark 
propagator S is known numerically. 

Note that in all bound state equations the real 
momentum-integration variables and the imaginary total 
on-shell momentum of the bound state are combined in the 
(anti-)quark momenta q± and q\, q 2 , q 3 to complex four- 
vectors. As a consequence, the arguments of the dressing 
functions in the quark propagator become complex. In the 
meson BSE, for example, the dressing functions are needed 
in a parabolic region in the complex plane defined by q±. 
Over the past years reliable numerical approaches to this 
problem have been developed and the required computa- 
tions are well under control (for more details see [74 - 76 1 V 

In this way, with the specification of the truncation used, 
one both decides the structure of the interaction kernel and 
obtains the quark propagator consistent with this kernel. 

3. Numerical representation 

The first step towards a numerical representation of 
Eqs. ([1]), @, and © is the analysis of the Lorentz and 
Dirac structure of the respective amplitudes. This struc- 
ture is a result of the particular representation properties 
of the state under consideration under the Lorentz group, 
including the state's parity and spin. Therefore, the bound 
state amplitudes are decomposed into Lorentz-covariant 
parts Ti and Lorentz- invariant parts F l , respectively, read- 
ing 



N 



(4) 



where the number of terms N as well as the tensor struc- 
ture of the Ti and V depend on the quantum numbers of 
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the bound state and all arguments have been suppressed 
for simplicity. The T are usually referred to as covariants, 
whereas we call the F 1 the components of the amplitude 
r. The Ti represent a spin-determined basis for the bound 
state, and one is — to some extent — free to choose the 
details thereof. 

To be more concrete, we consider the case of mesons 
in more detail, where two spin-1/2 fermions (the quark 
and antiquark) are combined to a boson with total spin J. 
As a result one obtains a 4 x 4-matrix structure with the 
correct Lorentz-transformation properties for a particle of 
spin J, see e.g. [37j [ . Consequently one has, in addition 
to the total meson momentum P M and the relative mo- 
mentum g M between the constituents, another four vector 
7 M of Dirac matrices to construct the elements of the BSA. 
The various combinations of the momenta and 7 M correctly 
encode the angular momentum structures inside the me- 
son, i.e. the contributions of (anti-)quark spin and orbital 
angular momentum. 

A scalar meson BSA, for example contains all possible 
Lorentz-scalar combinations of the three four-vectors P, q, 
and 7 (the actual construction is given below for a pseu- 
doscalar meson in Eqs. ([T7|) - ([20|) ). For the moment we 
note that while T and the Tj in Eq. (0| in general de- 
pend on the three four-vectors P, q, and 7 as such (which 
we denote by semicolons between arguments of an expres- 
sion), the components, being Lorentz- and Dirac-scalars, 
can only depend on scalar products of the momenta in- 
volved, i.e. P 2 , q 2 , and q ■ P in the meson case. Thus, 
writing arguments explicitly, Eq. Q reads 

JV 

r( 7 ; q; P) = £ T t ( T , <?; P) F l (P 2 ,q 2 ,q ■ P) , (5) 
i=l 

where again the tensor structure of the T, (and correspond- 
ingly r) was not denoted explicitly, since it is irrelevant to 
the following argument. 

In a numerical study, it is mostly advantageous to use a 
basis that is orthonormal, meaning the covariants satisfy 



Tr (Ti -Tj) = 5, 



(6) 



which also defines a generalized scalar product on the space 
of 4 x 4 matrices (any occurring Lorentz indices are un- 
derstood to be summed over here). Note that, for a set 
of covariants which is neither orthogonal nor normalized, 
the following step is more involved and we detail it in 
|Appendix C| If one uses the decomposition Q , the bound 
state equations Jl]), @, and ^ can be rewritten as cou- 
pled integral equations of the components depending on 
the scalar products of the momenta via the corresponding 
projections on the basis T. 

More concretely, we consider the integrand in, e.g. the 
(in-)homogeneous meson BSE, 

K [2] ( T .k-q-P) S a ( r ,q;P) r [2] ( T ,q;P) S b ( r ,q;P) ■ (7) 

The amplitude Tt 2 ] (7; q; P) expanded in the chosen Dirac 
basis Tj(j; q; P) and the result is projected on ^(7; k; P). 



Doing so, one obtains a matrix structure in the space of 
covariants, and Eq. ([7]) can be written as a matrix- vector 
multiplication in this space involving the BSE kernel ma- 
trix K)(k\q;P): 

K)(k-q ] P)F\P 2 ,q 2 ,q-P) = 

Tr [Ti( r , k; P) K [2] ( T , k; q; P) S a ( T , q; P) 
X T 3 ( T ,q;P) S b ( T ,q;P)]Fi(P 2 ,q 2 ,q- P) , (8) 

where the sum over the repeated index j is implied. 

The index j of the components F^(P 2 ,q 2 ,q ■ P) can 
thus be viewed as a vector index, which has to be con- 
tracted with the corresponding index of the kernel matrix 
Kj(k;q;P). Note that this procedure, although exempli- 
fied here for the case of mesons, is completely general, i.e., 
it applies to baryons as well and is valid for any choice of 
the interaction kernel K. 

The next step is to make the dependence on the con- 
tinuous momentum variables P 2 , q 2 , and q- P numerically 
accessible. To achieve this, we apply the so-called Nystrom 
or quadrature method (cf. Chap. 4]), which amounts 
to replacing an integral by a sum over suitable quadrature 
weights and points and neglecting the error term. Ap- 
plying this method discretizes the integration variables, 
and consequently also the momentum dependence on the 
left hand side. The homogeneous and the inhomogeneous 
bound state equations can then be written as matrix equa- 
tions in the covariants and the discretized momenta and 
read 

in the homogeneous case, and 



F i,v = F, 



i,V 



(10) 



in the inhomogeneous case. The indices i, j label the com- 
ponents, the multi- indices V, Q stand for all discretized 
momentum variables (summation over repeated indices is 
implied). The matrix K = ^C'g is the same in both equa- 
tions, and subsumes the interaction kernel, the dressed 
propagators of the constituents, the Dirac- and Lorentz 
structure, as well as the discretized integrations. It is ap- 
plied to a vector F l ' V representing the homogeneous or 
inhomogeneous bound state amplitude. 

As an alternative to the Nystrom method, one can ex- 
pand the momentum dependence of the components into 
suitable sets of orthogonal functions, which can then be 
integrated. In this approach, the index V of the vector 
F 1 -' v contains the coefficients of the expansion rather than 
the values of the components at certain points in momen- 
tum space (for applications in the present context, see 

e.g. mM)- 

A partial application of this alternative is the use of a 
Chebyshev expansion of the dependence in an angle vari- 
able as described in |Appcndix A where one only keeps 
a finite number of Chebyshev moments in the representa- 
tion of the amplitude. This step has been widely used in 
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DSE studies of hadron spectra and properties, and the 
fidelity of the approximation investigated in detail, see 



e -g- |4J, ISO, |8l|- While for studies of hadron masses a 



few moments are sufficient, more are required in situa- 
tions where considerable changes of the frame of reference 
are needed, such as form factor calculations at large mo- 
mentum transfer [82| ; ultimately, in these situations the 
approximation needs to be abandoned 
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4. Solution methods 

4-1. Homogeneous equations: eigenvalue algorithms 

With the results of the preceding section, the homoge- 
neous bound state equation (BSE or Faddeev equation), 
given in Eq. (|9]) in index notation, can be written as 



F[h) = K • F [h] 



(11) 



using matrix-vector notation. As already mentioned in 
Sec. [U this equation is only valid at the on-shell points of 
the bound states in the respective channel, i.e. at certain 
values of the total momentum squared P 2 = —M 2 , where 
n = 0,1,2,... numbers the ground- and all excited states 
in the channel. To find such a value of P 2 , one investigates 
the spectrum of K as a function of P 2 , since Eq. (TTT|) cor- 
responds to an eigenvalue equation (with the dependence 
on P 2 made explicit) 



X(P 2 )F [h] (P 2 ) = K(P 2 )-F [h] (P 2 ) , 



(12) 



where the eigenvalue A(P 2 ) = 1. In other words, to nu- 
merically approach a solution of the equation, a part of 
the result has to be already known, namely the values 
M 2 , or — more precisely — the mass of the state one is 
looking for. The way out is a self-consistency argument, 
where the eigenvalue spectrum is plotted as a function of 
P 2 and those points with X n (P 2 ) = 1 are identified: the 
largest eigenvalue determines the ground state, the smaller 
ones in succession the excitations of the system (see also 
Fig. U below) . Typically, one is interested in roughly up 
to five eigenvalues, since higher excitations are both not 
well-enough understood in theory and hard to access ex- 
perimentally. 

A great variety of algorithms is available to numerically 
tackle these kinds of problems, and the most commonly 
used is a simple iterative method. Similar to the other 
algorithms discussed in this section, it relies on the multi- 
plication of the matrix K on a vector and can successively 
be applied to find also excited states, by projecting on 
states already obtained, see e.g. [84|. This simple method, 



however, is not able to resolve pairs of complex conjugate 
eigenvalues, which may, for example, occur in the meson 
BSE [85]. In addition, the total number of required matrix- 
vector multiplications increases for every additional eigen- 
value, as demonstrated in Sec. [5j 

These difficulties are overcome by the use of more ad- 
vanced algorithms. For this purpose, we use the implicitly 



which is frequently ap- 
87| . An application of 



restarted Arnoldi factorization [86( 
plied in lattice QCD studies, e.g. 
this algorithm to bound state calculations is demonstrated 
in Sec. [5j where we use it to solve the pseudoscalar-meson 
BSE and compare the efficiency of both methods. 

4-2. Inhomogeneous equations: matrix inversion 

In the most compact notation, the inhomogeneous 
bound state equation can be written as 



F(P 2 ) = F (P 2 ) + K(P 2 ) • F(P 2 ) 



(13) 



where the matrix K(P 2 ) is identical to the one in Eq. 
(fTTj) . and the vector Fo is given by the decomposition of 
r according to Eq. (gj, T° = J2t T * F o together with the 
discretization of a possible momentum dependence. 

Again, the simplest method to treat this problem is a 
direct iteration. Mathematically, this corresponds to the 
representation of the solution by a von Neumann series (cf. 
[771 . Chap. 4]), which can be shown to converge as long as 
the norm of the operator K is smaller than one, ||K|| < 1. 
For matrices, this norm can be related to the largest eigen- 
value, such that for P 2 > —Mq, the iteration converges. 
When P 2 approaches the ground state position — Mq from 
above, the number of iterations necessarily grows, and no 
convergence is obtained if P 2 < —M§, as demonstrated in 
Sec. El 

However, a solution is possible for any P 2 if one rewrites 
Eq. (|SJ as 



F = (1 - K) • F 



(14) 



i.e., F is given by the inhomogeneous term Fq multiplied 
by the matrix inverse of (1 — K). F can then be com- 
puted by e.g. inverting the matrix exactly, which has been 
successfully used to resolve bound-state poles in the in- 
homogeneous amplitude, as shown in [66| in the case of 
mesons. On the downside, the direct inversion of a matrix 
is computationally expensive, and it is not straightforward 
to parallelize the procedure. 

A better approach is to view Eq. (IT4)) as a linear sys- 
tem whose solution is to be found. Equations like this are 
very common and several algorithms have been developed 
for their solution. In particular, if the matrix (1 — K) is 
big, Eq. ([T4l is a typical application for the so-called Con- 
jugate Gradient (CG) algorithms. Many types of these 
iterative Krylov-space methods are available. In the case 
of the bound-state equations considered here, the matrices 
involved are neither hermitian nor symmetric, such that a 
good choice is the well-known Bi-Conjugate-Gradients sta- 
bilized (BiCGstab) algorithm^ , which is widely used for 
example in lattice QCD (cf. [13, Chap. 6.2], where also the 
algorithm is presented in detail) . 

5. Application: numerical solution of the meson 
BSE 

As an illustration, we apply the algorithms discussed 
above to solve the homogeneous and inhomogeneous 
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pseudoscalar-meson BSEs and compare their efficiency 
in terms of the number of matrix-vector multiplications 
needed to achieve a specified accuracy. For bigger prob- 
lems like baryons in a three-quark setup, the kernel matrix 
typically does not fit into memory, and thus has to be re- 
computed on the fly in each iteration- or matrix-vector 
multiplication step. In this case, one matrix-vector mul- 
tiplication is rather time consuming and it is desirable to 
keep the number of necessary multiplications as small as 
possible. 

For our test case here, however, we study the 
pseudoscalar-meson BSE, where the kernel matrix is small, 
but one can still investigate the questions at hand. We 
employ the rainbow-ladder truncation, i.e. the rainbow ap- 
proximation in the quark propagator DSE together with 
a ladder truncation of the corresponding quark-antiquark 
BSE. We define the difference in relative momenta k — q=: 
£ and the transverse projector with respect to the momen- 
tum q as T^ u (q) := (i^ — "^r 1 )- I R this truncation, the 

kernel K[ 2 ]{t> fc ; 95 p ) 01 tne BSEs, Eqs. JTJ and © is then 
given by 



K [2] ( r ,k;q;P) 



-7m" 



(15) 



where the effective interaction as a function of the 



momentum-squared s, introduced in Ref. 45j, reads 



V(s) = D —se 



(16) 



The term Tuv{s) implements the perturbative run- 
ning coupling of QCD for large s, preserving the one- 
loop renormalization-group behavior of QCD. The Gaus- 
sian term models the enhancement in the intermediate- 
momentum regime necessary to produce a reasonable 
amount of dynamical chiral symmetry breaking. It con- 
tains the parameters of the model, D and w, describing the 
overall strength and momentum-space width (correspond- 
ing to an inverse effective range) of the interaction. The 
behavior of the effective interaction in the far infrared is 
expected to be of minor relevance to ground-state proper- 
ties (cf. [8{| and references therein) . For the present study, 
we make a common choice for the parameters, namely 
D = 0.93 GeV 2 and oj = 0.4 GeV (for full details on the 
truncation, the effective coup ling, or the effects of other 
parameter values, see e.g. |39i. 144 l45j). 

5.1. Kernel setup 

The above definitions together with the dressed quark 
propagators computed already completely specify the in- 
gredients of the BSE. We investigate light quarks in anal- 
ogy to [39j |. where the current-quark mass in an isospin- 
symmetric setup was adjusted to fit the bound-state mass 
of the p meson. The details of the discretization of the 
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Figure 4: The five largest eigenvalues of the homogeneous BSE plot- 
ted over \/— P 2 which corresponds to the bound state mass M where 
A = 1 indicated by the horizontal dashed line. The ground-state 
(leftmost intersection) solution vector has positive C-parity (pion), 
the second has negative (exotic) and the third again has positive 
C-parity (excited pion). 



kernel matrix proceed as follows: The orthonormal pseu- 
doscalar covariants, constructed to satisfy Eq. ((6]), read 



75 

2 ' 



T, = 



Ta = 



I^P 2 

75((7-< z)-^^) 

^7s((7-g)(7--P)-(7--P)(7-g)) 
2 V /P 2 g 2 - (P ■ q) 2 



(17) 
(18) 

(19) 
(20) 



Choosing the rest frame of the quark-anti-quark system, 
and applying the parametrization and discretization as de- 
scribed in |Appcndix A[ the kernel matrix Eq. @ in our 
setup becomes 

4 r 1 v (f) 

- p»W Jjy -Y^d) 

x Tr [T i ( r ,k;P)^S(q + )T j ( r ,q;P)S(q-)^} , (21) 

where w[q\ f], w[z m ] denote the quadrature weights and the 
replacements k 2 — > fc 2 , z% — ?> z s , q 2 — > kf , z — > z m have 
been made in all occurring momenta to implement the 
discretization. Therefore, the indices i;j label the compo- 
nents and r, s; m the momentum space points. For the 
following calculations, we use N q = 32 and N z — 24, such 
that K has the dimensions (32,24,4) x (32,24,4). 



5.2. Homogeneous BSE 

To solve the homogeneous BSE, we use both the MPI 
based version of the ARPACK library (an implementation 
of the implicitly restarted Arnoldi factorization) and the 
simple iteration. Fig. @] shows the largest five eigenvalues 
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Figure 5: The number of matrix- vector multiplications needed for 
convergence of the simple iteration (circles) and the Arnoldi factor- 
ization (squares), plotted against the number of eigenvalues com- 
puted at a (typical) fixed value of P 2 = —M 2 . 
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Figure 6: Component F±(P 2 ,0, 0) of the inhomogeneous pseu- 
doscalar amplitude calculated using BiCGstab vs. the square of the 
total momentum P 2 . The vertical lines mark the pole positions, cor- 
responding to the pion ground- and first excited state (J PC = + ~ ). 



of K for our example. Bound-state masses can be identi- 
fied by the positions at which an eigenvalue curve crosses 
one (dashed line in the figure) and from left to right cor- 
respond to the ground- and first-excited, second-excited, 
etc., states. Note that in our approach we do not restrict 
the symmetry of the amplitudes (eigenvectors) with re- 
spect to the the angular variable z, such that we obtain ho- 
mogeneous solutions of both positive and negative charge- 
conjugation parity (C-parity) for equal-mass constituents 
and a choice of 1/2 for the momentum-partitioning pa- 
rameters r]±, as indicated above (for more details on the 
definition and calculation of C, see |Appendix B[ ) . In the 



pseudoscalar case, a negative C-parity is considered exotic, 
since it is not available for a qq state in quantum mechan- 
ics. However, it appears naturally in a quark-antiquark 
BSE setup, where our main interest here comes from a sys- 
tematic point of view. A more general discussion of states 
with exotic C-parity in this formalism and their possible 
interpretations can be found, e.g., in [37, 

To compare the efficiency of the two algorithms, Arnoldi 
factorization versus iteration, we compute the one to four 
largest eigenvalues of K and compare the convergence in 
terms of the number of iterations needed to obtain an ab- 
solute accuracy of the eigenvector of e = 10~ 8 , at a typical 
value of P 2 = -Af 2 = 0.0527 GeV 2 . The results, given 
in Fig. [5j show that for the first eigenvalue the Arnoldi 
factorization needs only half as many matrix- vector multi- 
plications as the simple iteration. With increasing number 
of eigenvalues, the use of this advanced algorithm becomes 
even more advantageous. 

Another interesting observation from Fig. [5] is that the 
Arnoldi factorization was more efficient for three eigenval- 
ues than when only two were requested. This is most likely 
due to a "clustering" of eigenvalues number two and three 
for the algorithm, an effect which appears for eigenvalues 



close together and is also related to the eigenvectors. In 
this particular case, eigenvectors two and tree have oppo- 
site C-parity or z-symmetry, which may make them more 
easily distinguishable for the algorithm and more easy to 
obtain as a result. The ARPACK library is very efficient 
at evaluating all eigenvalues in such a cluster, while con- 
vergence is slower, if one asks for only one or a few of 
the ei genv alues in the cluster (see also the ARPACK users 
guide [92|). 

5.3. Inhomogeneous BSE 

We apply the direct iteration (summation of the von 
Neumann series) and the inversion using the BiCGstab al- 
gorithm to solve the inhomogeneous BSE ((TJ), in the setup 
described above for pseudoscalar quantum numbers. 

In the inhomogeneous case not only the structure of the 
amplitude determines the quantum numbers of the solu- 
tion but also that of the inhomogeneous term Tq. Follow- 
ing [gij], a possible choice for pseudoscalars is 



r = Z475 



(22) 



where Z4 is a renormalization constant obtained from the 
gap equation (cf. [44]). With this choice (pseudoscalar, 
positive C-parity), no poles corresponding to negative C- 
parity appear in the solution, as can be seen from Fig. [SJ 
The curve shown in this figure has been obtained with 
the BiCGstab algorithm, because as described in Sec. 14.21 
the direct iteration fails to converge if P 2 < —Mq. This is 
demonstrated in Fig. [3 where the number of matrix- vector 
multiplications needed for convergence is plotted against 
P 2 for both methods. 

It is clear that the number of matrix-vector multipli- 
cations needed for the direct iteration diverges as P 2 ap- 
proaches — Mq (note that Fig.[7]uses a logarithmic scale on 
the vertical axis). The inversion with BiCGstab, however, 
converges for all P 2 with nearly the same speed, needing 
approximately 10 matrix-vector multiplications. 
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BiCGstab algorithm (squares), plotted on a logarithmic scale against 
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line indicates the position of the ground state P 2 = —Mq of the 
system. Note that the straightforward iteration does not converge 
for P 2 < -M 2 . 



6. Conclusion and Outlook 



We have investigated aspects and benefits of different 
numerical solution methods for the pseudoscalar meson 
BSE in a rainbow-ladder truncation of the DSE approach 
to QCD. In a realistic model, our comparison was aimed 
at a small-scale test and subsequent identification of ef- 
ficient algorithms for the numerical approach to bound- 
state problems in this covariant continuum approach to a 
quantum field theory. Since the algorithms used and ad- 
vocated are applicable in a very general way to the matrix 
representations of the integral kernels appearing in such 
bound-state studies, the strategies proposed here are valu- 
able for similar studies of bound states in more involved 
truncations than rainbow-ladder on one hand. 

On the other hand, any bound-state problem involving 
more than two constituents, starting with but not limited 
to baryons in a three-quark setup, needs efficient meth- 
ods to perform sophisticated numerical studies, given the 
computational resources presently available. In general 
one will find that in both examples given here one simply 
has to deal with a larger kernel matrix and thus efficiency 
of the algorithms used is of the essence. 

One particularly interesting observation with regard to 
present-day covariant three-quark studies is that comput- 
ing bound-state masses is more efficient using the inhomo- 
geneous equation than the homogeneous one. In addition, 
the kernel matrices often need to be constructed on the 
fly due to limited system memory. With these two aspects 
combined, such a problem appears to be an ideal applica- 
tion for the field of GPU computing. 



Appendix A. Aspects of four-dimensional mo- 
mentum integration 

We use 4-dim. spherical coordinates, such that the mo- 
mentum integration is written as 



^ d(q 2 )^ J dzVT^~z 2 / ,/.,/ I ,/,,. 



(A.1) 

In the meson BSE, there are three relevant momenta: P 
(total momentum), k (relative momentum), and q (loop 
momentum) . Subsequently, we choose P to be in the rest- 
frame of the bound state, 

P= (0,0,0,VP2) . (A.2) 
The other momenta are chosen accordingly and read 



k 2 ( 0,0, x /l-z 2 ,z k 



(A.3) 



and 



q = 



(o, Vi-zVi-y 2 ,\/i-^,z) • (A.4) 



In this parametrization, the integration J d<p is trivial. 

The components of the amplitude on the left hand side 
of the BSEs Eqs. ([1) and @ depend on the scalar products 
P ■ k = ZkVk^VP 2 , P 2 , and k 2 . Inside the integral, on 
the right hand side, the scalar products become P ■ q = 
zy^fVP 2 , P 2 , and q 2 . Thus, the BSE kernel matrix K 
induces the following mapping on the momentum variables 



(q 2 ,z) h-> (k 2 ,z k ) 



(A.5) 



such that the integration J dy does not add a dimension 
K, although it is not trivial. 

After choosing a parametrization of the momentum in- 
tegration, the next step is to discretize the momentum 
dependence. In this work, we straightforwardly apply the 
quadrature method, and replace 



d{q 2 ) 



Na N z 



dz\f\ 



£ £ w[qf]w[z m ] , 



1=1 m=l 



where w[qf], 



9i 



(A.6) 

2 2 n ] denote the quadrature weights and 
the corresponding nodes. The factors of q 2 jl and 



VI — z 2 have been absorbed in the weights. 

Note that, especially for the integration over z, it is ad- 
vantageous to use a quadrature rule whose weights include 
VT — z 2 by construction, e.g., the Gauss-Chebyshev type 
2 rule. 

An alternative, advantageous and widely used in the 
calculations of hadron spectra is to apply the quadrature 
method discussed above only to J d(q 2 ) and to resolve the 
z-dependence and by an expansion in Chebyshev polyno- 
mials of the second kind. The components are then written 
as 



F\P 2 ,q 2 ,z) 



M 

£ 

m— 1 



l Fi(q 2 ,P 2 )U m (z) 



(A.7) 
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with m Fi(q 2 , P 2 ) the so-called Chebyshev moments, which 
retain only the functional dependence on k 2 and P 2 . The 
number of terms M taken into account is finite in practice, 
but infinite in principle. The Chebyshev polynomials of 
the second kind U m (z) satisfy the orthogonality relation 

2 J dzVl - z 2 U m {z) U n {z) = 6 mn . (A.8) 

To obtain a matrix structure not only in the covariants, 
but also in the Chebyshev moments, the above expansion 
is inserted in Eq. (j5J), and is then projected on one moment 
by use of Eq. (IA.8I) . This finally leads to 

J K*£(k;q;P) n Fj( q 2 ,P 2 ) = 

J J rfz fc f/ m (z fe )Tr[r i (7;fc;P)K [2] ( 7 ;fc;( Z ;P) x 

s a ( r , gjP^-fygP) s b ( T ,q;P)] u n (z)) n F 3 { q 2 ,p 2 ). 

(A.9) 

Again, the sum over the repeated indices j, n is implied. 

Appendix B. (Generalized) C-parity 

The action of the C-parity transformation on the meson 
BSA is defined via 

T(P;k) = (CT(P;-k)C- 1 ) t , (B.f) 

where t denotes the matrix-transpose, and the charge- 
conjugation matrix C = 7274- If the C-parity is a good 
quantum number of the system, T(P; k) is an eigenstate 
of the C-parity operation given above, with eigenvalues 
Ac = ±l, 

(Cr(P; -fcJCT 1 )* = A c r(P; k) . (B.2) 

The amplitudes occurring in Eq. (|B.2|) can be decomposed 
into covariants and components according to Eq. (j5j) , such 
that in the rest frame of the bound state with the notation 
introduced in Eqs. (|A.2I) - (|A.4I) . 

N 

J2(C T,(P; -k) C- 1 ) 1 F\P 2 , k 2 , -z k ) 
i=l 

N 

= \ c Y,Ti{P;k)F\P 2 ,k 2 ,z k ) ■ (B.3) 

i=l 

Projecting this equation on one covariant using the orthog- 
onality relation ([6]), we obtain 

Fi(P 2 ,k 2 ,z k ) := 

N 

Y J ^[T ] {P;k){CT l {P--k)C- 1 ) t ] F\P 2 ,k 2 ,-z k ) 

i=l 

= \ c Fi(P 2 ,k 2 ,z k ) . (B.4) 



The next step is to discretize the momenta using the 
quadrature method, such that the functional dependence 
of the component vectors on the momentum variables can 
be represented in index notation, 

Fi(P 2 ,k 2 ,z k ) => F^ m , (B.5) 
P^P 2 ,*; 2 ,^) j?M,m, ( B6 ) 

Now, the amplitudes are given as complex vectors, such 
that the canonical scalar product on C" can be used to 
solve for Ac, 

Ac = — — tft- — f- , (B.7) 

where * denotes complex conjugation, and repeated indices 
are summed over. 

For states with definite C-parity Ac = ±1- As can be 
seen from the above equations, the C-parity is determined 
by the (anti-)symmetry of the amplitudes with respect to 
z k . For states which are not eigenstates of the C operation, 
Eq. (|B.7|) can be used to define a " generalized C-parity" . 
It lies between —1 and 1, and its deviation from these 
values indicates the asymmetry of the state caused, e.g., 
by mass difference of the constituents. 

It is interesting to note that the C-parity as a symme- 
try property of the eigenvectors of K is constant over the 
whole range of P 2 , even if the on-shell condition for this 
state is not fulfilled, i.e. P 2 ^ —M 2 . 

Appendix C. Using a non-orthogonal Dirac basis 

Consider the inhomogeneous BSE written in the general 
form 

/( 7 ;fc;P)r( 7 ;A : ;P) = Z(7;fc;P) 

- f K^r.k-q-P^q-P^r.k-q-.P) , (C.l) 

where the dependence of every term on all variables includ- 
ing 7 matrices is given explicitly. The ; between variables 
again denotes a dependence on complete four- vectors. K\ 
and Ki represent generalized formal kernel pieces, Z a 
general driving term, and / an arbitrary function of its ar- 
guments. To transform this equation into a set of coupled 
integral equations for components and then use Chebyshev 
moments, which are described in |Appcndix A| we write 
the BSA as the sum over its covariants Tj and Lorentz- as 
well as Dirac-scalar components F % and the latter as sums 
over Chebyshev polynomials and moments 

M N 

r( 7 ; k; P) = T ^ *; P) l F 3 (k 2 , P 2 )Ui(z k ) , (C.2) 

1=1 J=l 

where M is the number of Chebyshev polynomials taken 
into account and N is the number of covariants in the 
BSA. The Chebyshev polynomial of the second kind Ui(z k ) 
depends on z k := k ■ Pj \Jk 2 P 2 . Now we apply T n (j; k; P) 
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on Eq. (|C.1|) from the left and take the Dirac trace. The 
result is 

M N 

J2 A «^ fc2 ' p2 ' z kYFj(k 2 , P 2 )Ui(z k ) = 

i=l j=l 

M N 

Z n (k\P\z k )-Y,Yl / Tr [T„( 7 ; k; P) 

1=1 m=l J 1 

x Ki(rM*P)T m {a\*\P) l F m { q \p 2 ) 

xC/ i (z 9 )^ 2 ( 7 ;fc; g ;P)] , (C.3) 

where 

A nj {k\ P 2 lZk ) := Tr [T n ( T , k; P)/( 7 ; k; P)T,( 7 ; k; P)] 
Z n (k 2 ,P 2 ,z k ) :=Tt[T n {r,k;P)Z{r,k;P)\ (C.4) 

The next step is to invert the matrix A n j for each set 
of coordinates (k 2 ,P 2 ,z k ), and apply its inverse to the 
equation, i.e., J2n=i -^rn fr° m the left: 

M N 

]T l F r (k 2 , P 2 )Ui(z k ) = A~n(k 2 ,P 2 , z k )Z n (k 2 , P 2 ,z k ) 

i=l n=l 
N M N 

"EEE / A- 1 (fc 2 ,P 2 ,z fc )Tr[T„(7;fc;P)K 1 ( 7 ;fc;g;P) 

n=l 1=1 m=l J 1 

x T m (7;< Z ;P) i F m (g 2 ,P 2 )^(z g )^ 2 ( 7 ;fc;g;P)] (C.5) 

The last step is the projection with the help of Cheby- 
shev polynomials via — J dz k y/ 1 — z\ Uj(z k ) from the left 
(cf. Eq. (|A.8[) 1 and one obtains 

M N 

iF r (k 2 ,p 2 )=v Z (j, r ,k 2 )-Y: e - / / A 7 ^ 

i=l n,m=l 77 

x Uj(z k )A^(k 2 , P 2 , z fc )Tr [T n (7j fc; P)^i( 7 ; fc; «; P) 
xT m ( 7 ;g;P) ; P m (g 2 , P 2 )^(^)^ 2 ( 7 ; fc; <z; P)] (C.6) 
with the driving term given by 

2 N r i 

V z (j, r, k 2 ) := - Y / ^V 1 ~ z fc 

n=l J 

x ^(fc^J* ZfcjZ^.P 2 ,**) . (C.7) 

This procedure does not require the set of covariants to 
be orthogonal, it is completely general, also with respect 
to the kernel, the driving term, and possible terms multi- 
plying the amplitude on the left-hand side of the BSE. All 
terms and projections are included correctly via the matrix 
A. In the case considered in the present work, / = 1 and 
the driving term has the standard form for pseudoscalar 
mesons, Eq. ([22]). 
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